home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT8.CPP < prev    next >
C/C++ Source or Header  |  1995-01-16  |  9KB  |  364 lines
  1. //$$ newmat8.cpp         Advanced LU transform, scalar functions
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmatap.h"
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,8); ++ExeCount; }
  12.  
  13. #define REPORT {}
  14.  
  15.  
  16. /************************** LU transformation ****************************/
  17.  
  18. void CroutMatrix::ludcmp()
  19. // LU decomposition - from numerical recipes in C
  20. {
  21.    REPORT
  22.    Tracer trace("Crout(ludcmp)");
  23.    int i,j;
  24.  
  25.    Real* vv=new Real [nrows]; MatrixErrorNoSpace(vv);
  26.    MONITOR_REAL_NEW("Make  (CroutMat)",nrows,vv)
  27.    Real* a;
  28.  
  29.    a=store;
  30.    for (i=0;i<nrows;i++)
  31.    {
  32.       Real big=0.0;
  33.       j=nrows; while (j--) { Real temp=fabs(*a++); if (temp > big) big=temp; }
  34.       if (big == 0.0)
  35.       {
  36.          MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  37. #ifdef Version21
  38.          sing = TRUE; delete [] vv; return;
  39. #else
  40.          sing = TRUE; delete [nrows] vv; return;
  41. #endif
  42.       }
  43.       vv[i]=1.0/big;
  44.    }
  45.  
  46.    Real* aj=store;
  47.    for (j=0;j<nrows;j++)
  48.    {
  49.       Real* ai=store;
  50.       for (i=0;i<j;i++)
  51.       {
  52.          Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  53.          int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  54.          *ajx = sum; ai += nrows;
  55.       }
  56.  
  57.       Real big=0.0; int imax;
  58.       for (i=j;i<nrows;i++)
  59.       {
  60.          Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  61.          int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  62.          *aix = sum; ai += nrows;
  63.          Real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
  64.       }
  65.  
  66.       if (j != imax)
  67.       {
  68.          Real* amax=store+imax*nrows; Real* ajx=store+j*nrows;
  69.          int k=nrows; while(k--) { Real dum=*amax; *amax++=*ajx; *ajx++=dum; }
  70.          d=!d; vv[imax]=vv[j];
  71.       }
  72.  
  73.       indx[j]=imax; ai=aj+j*nrows;
  74.       if (*ai == 0.0)
  75.       {
  76.          MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  77. #ifdef Version21
  78.          sing = TRUE; delete [] vv; return;
  79. #else
  80.          sing = TRUE; delete [nrows] vv; return;
  81. #endif
  82.       }
  83.       Real dum=1.0/(*ai);
  84.       i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
  85.  
  86.       aj++;
  87.    }
  88.    MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  89. #ifdef Version21
  90.    delete [] vv;
  91. #else
  92.    delete [nrows] vv;
  93. #endif
  94. }
  95.  
  96. void CroutMatrix::lubksb(Real* B, int mini)
  97. {
  98.    REPORT
  99.    Tracer trace("Crout(lubksb)");
  100.    if (sing) Throw(SingularException(*this));   
  101.    int i,j; int ii=-1; Real* ai=store;
  102.  
  103.    for (i=0;i<nrows;i++)
  104.    {
  105.       int ip=indx[i]; Real sum=B[ip]; B[ip]=B[i];
  106.       if (ii>=0)
  107.       {
  108.          Real* aij=ai+ii; Real* bj=B+ii; j=i-ii;
  109.          while (j--) { sum -= *aij++ * *bj; bj++; }
  110.       }
  111.       else if (sum) ii=i;
  112.       B[i]=sum; ai += nrows;
  113.    }
  114.  
  115.    for (i=nrows-1;i>=mini;i--)
  116.    {
  117.       Real* bj=B+i; ai -= nrows; Real* ajx=ai+i; Real sum=*bj; bj++;
  118.       j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
  119.       B[i]=sum/(*(ai+i));
  120.    }
  121. }
  122.  
  123.  
  124. /****************************** scalar functions ****************************/
  125.  
  126. inline Real square(Real x) { return x*x; }
  127.  
  128. Real GeneralMatrix::SumSquare() const
  129. {
  130.    REPORT
  131.    Real sum = 0.0; int i = storage; Real* s = store;
  132.    while (i--) sum += square(*s++);
  133.    ((GeneralMatrix&)*this).tDelete(); return sum;
  134. }
  135.  
  136. Real GeneralMatrix::SumAbsoluteValue() const
  137. {
  138.    REPORT
  139.    Real sum = 0.0; int i = storage; Real* s = store;
  140.    while (i--) sum += fabs(*s++);
  141.    ((GeneralMatrix&)*this).tDelete(); return sum;
  142. }
  143.  
  144. Real GeneralMatrix::Sum() const
  145. {
  146.    REPORT
  147.    Real sum = 0.0; int i = storage; Real* s = store;
  148.    while (i--) sum += *s++;
  149.    ((GeneralMatrix&)*this).tDelete(); return sum;
  150. }
  151.  
  152. Real GeneralMatrix::MaximumAbsoluteValue() const
  153. {
  154.    REPORT
  155.    Real maxval = 0.0; int i = storage; Real* s = store;
  156.    while (i--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
  157.    ((GeneralMatrix&)*this).tDelete(); return maxval;
  158. }
  159.  
  160. Real SymmetricMatrix::SumSquare() const
  161. {
  162.    REPORT
  163.    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  164.    for (int i = 0; i<nr; i++)
  165.    {
  166.       int j = i;
  167.       while (j--) sum2 += square(*s++);
  168.       sum1 += square(*s++);
  169.    }
  170.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  171. }
  172.  
  173. Real SymmetricMatrix::SumAbsoluteValue() const
  174. {
  175.    REPORT
  176.    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  177.    for (int i = 0; i<nr; i++)
  178.    {
  179.       int j = i;
  180.       while (j--) sum2 += fabs(*s++);
  181.       sum1 += fabs(*s++);
  182.    }
  183.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  184. }
  185.  
  186. Real SymmetricMatrix::Sum() const
  187. {
  188.    REPORT
  189.    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  190.    for (int i = 0; i<nr; i++)
  191.    {
  192.       int j = i;
  193.       while (j--) sum2 += *s++;
  194.       sum1 += *s++;
  195.    }
  196.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  197. }
  198.  
  199. Real BaseMatrix::SumSquare() const
  200. {
  201.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  202.    Real s = gm->SumSquare(); return s;
  203. }
  204.  
  205. Real BaseMatrix::SumAbsoluteValue() const
  206. {
  207.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  208.    Real s = gm->SumAbsoluteValue(); return s;
  209. }
  210.  
  211. Real BaseMatrix::Sum() const
  212. {
  213.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  214.    Real s = gm->Sum(); return s;
  215. }
  216.  
  217. Real BaseMatrix::MaximumAbsoluteValue() const
  218. {
  219.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  220.    Real s = gm->MaximumAbsoluteValue(); return s;
  221. }
  222.  
  223. Real Matrix::Trace() const
  224. {
  225.    REPORT
  226.    Tracer trace("Trace");
  227.    int i = nrows; int d = i+1;
  228.    if (i != ncols) Throw(NotSquareException(*this));
  229.    Real sum = 0.0; Real* s = store;
  230.    while (i--) { sum += *s; s += d; }
  231.    ((GeneralMatrix&)*this).tDelete(); return sum;
  232. }
  233.  
  234. Real DiagonalMatrix::Trace() const
  235. {
  236.    REPORT
  237.    int i = nrows; Real sum = 0.0; Real* s = store;
  238.    while (i--) sum += *s++;
  239.    ((GeneralMatrix&)*this).tDelete(); return sum;
  240. }
  241.  
  242. Real SymmetricMatrix::Trace() const
  243. {
  244.    REPORT
  245.    int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  246.    while (i--) { sum += *s; s += j++; }
  247.    ((GeneralMatrix&)*this).tDelete(); return sum;
  248. }
  249.  
  250. Real LowerTriangularMatrix::Trace() const
  251. {
  252.    REPORT
  253.    int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  254.    while (i--) { sum += *s; s += j++; }
  255.    ((GeneralMatrix&)*this).tDelete(); return sum;
  256. }
  257.  
  258. Real UpperTriangularMatrix::Trace() const
  259. {
  260.    REPORT
  261.    int i = nrows; Real sum = 0.0; Real* s = store;
  262.    while (i) { sum += *s; s += i--; }
  263.    ((GeneralMatrix&)*this).tDelete(); return sum;
  264. }
  265.  
  266. Real BandMatrix::Trace() const
  267. {
  268.    REPORT
  269.    int i = nrows; int w = lower+upper+1;
  270.    Real sum = 0.0; Real* s = store+lower;
  271.    while (i--) { sum += *s; s += w; }
  272.    ((GeneralMatrix&)*this).tDelete(); return sum;
  273. }
  274.  
  275. Real SymmetricBandMatrix::Trace() const
  276. {
  277.    REPORT
  278.    int i = nrows; int w = lower+1;
  279.    Real sum = 0.0; Real* s = store+lower;
  280.    while (i--) { sum += *s; s += w; }
  281.    ((GeneralMatrix&)*this).tDelete(); return sum;
  282. }
  283.  
  284. Real BaseMatrix::Trace() const
  285. {
  286.    REPORT
  287.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(MatrixType::Dg);
  288.    Real sum = gm->Trace(); return sum;
  289. }
  290.  
  291. void LogAndSign::operator*=(Real x)
  292. {
  293.    if (x > 0.0) { log_value += log(x); }
  294.    else if (x < 0.0) { log_value += log(-x); sign = -sign; }
  295.    else sign = 0;
  296. }
  297.  
  298. Real LogAndSign::Value() const { return sign * exp(log_value); }
  299.  
  300. LogAndSign::LogAndSign(Real f)
  301. {
  302.    if (f == 0.0) { log_value = 0.0; sign = 0; return; }
  303.    else if (f < 0.0) { sign = -1; f = -f; }
  304.    else sign = 1;
  305.    log_value = log(f);
  306. }
  307.  
  308. LogAndSign DiagonalMatrix::LogDeterminant() const
  309. {
  310.    REPORT
  311.    int i = nrows; LogAndSign sum; Real* s = store;
  312.    while (i--) sum *= *s++;
  313.    ((GeneralMatrix&)*this).tDelete(); return sum;
  314. }
  315.  
  316. LogAndSign LowerTriangularMatrix::LogDeterminant() const
  317. {
  318.    REPORT
  319.    int i = nrows; LogAndSign sum; Real* s = store; int j = 2;
  320.    while (i--) { sum *= *s; s += j++; }
  321.    ((GeneralMatrix&)*this).tDelete(); return sum;
  322. }
  323.  
  324. LogAndSign UpperTriangularMatrix::LogDeterminant() const
  325. {
  326.    REPORT
  327.    int i = nrows; LogAndSign sum; Real* s = store;
  328.    while (i) { sum *= *s; s += i--; }
  329.    ((GeneralMatrix&)*this).tDelete(); return sum;
  330. }
  331.  
  332. LogAndSign BaseMatrix::LogDeterminant() const
  333. {
  334.    REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  335.    LogAndSign sum = gm->LogDeterminant(); return sum;
  336. }
  337.  
  338. LogAndSign GeneralMatrix::LogDeterminant() const
  339. {
  340.    REPORT
  341.    Tracer tr("Determinant");
  342.    if (nrows != ncols) Throw(NotSquareException(*this));
  343.    CroutMatrix C(*this); return C.LogDeterminant();
  344. }
  345.  
  346. LogAndSign CroutMatrix::LogDeterminant() const
  347. {
  348.    REPORT
  349.    if (sing) return 0.0;
  350.    int i = nrows; int dd = i+1; LogAndSign sum; Real* s = store;
  351.    while (i--) { sum *= *s; s += dd; }
  352.    if (!d) sum.ChangeSign(); return sum;
  353.  
  354. }
  355.  
  356. LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
  357. : gm( ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver() )
  358. {
  359.    if (gm==&bm) { REPORT  gm = gm->Image(); } 
  360.    // want a copy if  *gm is actually bm
  361.    else { REPORT  gm->Protect(); }
  362. }
  363.  
  364.